function [result, ...
          T, Tcoord,n_U_sets, n_vertices, C,idU1, idU2, idU3, idU4, s1,s2,s3,C_sign, sU, Aeq, beq]=...
          OneSide(param_MOSEK, g,id_Ugridreduced, s_id_gridreduced, sU_id_gridreduced,...
                  P_col,...
                  Tcoord_zeros,...
                  Ugridreduced, n_U_sets, indices_pairs_extended, Tcoord, C, C_sign,...
                  cr_obsequiv, cr_one, cr_zeros,fill_obsequiv, fill_one, fill_zeros, numeqconstraints, RHS_one, RHS_zeros,RHS_obsequiv,n_vertices,infty_pos,...
                  cr_symmetry, fill_symmetry, RHS_symmetry, S1, S2, S3,...
                  ind)      
% Set indices to keep track of equal elements     
idU1=id_Ugridreduced{1}(g,:);
idU2=id_Ugridreduced{2}(g,:); 
idU3=id_Ugridreduced{3}(g,:); 
idU4=id_Ugridreduced{4}(g,:); %for each of the U-sets
% Relevant sizes of the vector of unknowns
s1=s_id_gridreduced(g,1);
s2=s_id_gridreduced(g,2);
s3=s_id_gridreduced(g,3);
%s4=s_id_gridreduced(g,4);
sU=sU_id_gridreduced(g); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Equalities %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Observational equivalence
% Coordinate columns
cc_obsequiv=[changecoord_4(s1,s2,s3, idU1(P_col{1}(:,1)), idU2(P_col{1}(:,2)), idU3(P_col{1}(:,3)), idU4(P_col{1}(:,4))) ... %Y=1
             changecoord_4(s1,s2,s3, idU1(P_col{2}(:,1)), idU2(P_col{2}(:,2)), idU3(P_col{2}(:,3)), idU4(P_col{2}(:,4))) ... %Y=2
             changecoord_4(s1,s2,s3, idU1(P_col{3}(:,1)), idU2(P_col{3}(:,2)), idU3(P_col{3}(:,3)), idU4(P_col{3}(:,4))) ... %Y=3
             changecoord_4(s1,s2,s3, idU1(P_col{4}(:,1)), idU2(P_col{4}(:,2)), idU3(P_col{4}(:,3)), idU4(P_col{4}(:,4))) ... %Y=4
             changecoord_4(s1,s2,s3, idU1(P_col{5}(:,1)), idU2(P_col{5}(:,2)), idU3(P_col{5}(:,3)), idU4(P_col{5}(:,4)))];   %Y=0

%% Integrates to one
cc_one=changecoord_4(s1,s2,s3, idU1(infty_pos), idU2(infty_pos), idU3(infty_pos), idU4(infty_pos));

%% Zeros
cc_zeros=changecoord_4(s1,s2,s3, idU1(Tcoord_zeros(:,1)), idU2(Tcoord_zeros(:,2)), idU3(Tcoord_zeros(:,3)), idU4(Tcoord_zeros(:,4))); %1xn_zeros

%% Symmetry
cc_symmetry=[reshape([changecoord_4(s1, s2, s3, idU1(S1(:,1)), idU2(S1(:,2)), idU3(S1(:,3)), idU4(S1(:,4))).' ... 
                      changecoord_4(s1, s2, s3, idU1(S2(:,1)), idU2(S2(:,2)), idU3(S2(:,3)), idU4(S2(:,4))).'].', 1, (n_U_sets+1)*n_U_sets*2) ...
             ...         
             changecoord_4(s1,s2, s3, idU1(S3(:,1)), idU2(S3(:,2)), idU3(S3(:,3)), idU4(S3(:,4)))]; 
         
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Inequalities %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Increasingness 
%Repeat the way we constructed Tcoord but now using the parameter values under consideration
vectors_augmented = cellfun(@(x) {x(g,:)}, Ugridreduced); %cell 1xn_U_sets
%For j=1,...,n_U_sets, it picks the g-th row of Ugridreduced{j}
T_temp = cell(1,n_U_sets);
[T_temp{:}] = ndgrid(vectors_augmented{:}); 
T_temp = cat(n_U_sets+1, T_temp{:});
T_augmented = reshape(T_temp,[],n_U_sets); %matrix of dimension as Tcoord reporting the elements of all the possible n_U_sets-tuples from the U-sets
T = T_augmented(ind, :); %reorder so as to mimic Restrictions1

D=[T(indices_pairs_extended(:,1),:) T(indices_pairs_extended(:,2),:) ...
   Tcoord(indices_pairs_extended(:,1),:) Tcoord(indices_pairs_extended(:,2),:)]; %n_U_sets-tuples n_U_sets-tuples coordinates coordinates


%Save in D1 the rows of D where the first n_U_sets-tuple is <= than the second
%Only the columns reporting column-coordinates are relevant
D1=D(sum(D(:,1:n_U_sets)<=D(:,n_U_sets+1:2*n_U_sets),2)==n_U_sets, 2*n_U_sets+1:end); %sd1x(2*n_U_sets)
sd1=size(D1,1);

%Save in D2 the rows of D where the first n_U_sets-tuple is >= than the second
%Only the columns 5:end (reporting column-coordinates) are relevant
D2=D(sum(D(:,1:n_U_sets)>=D(:,n_U_sets+1:2*n_U_sets),2)==n_U_sets, 2*n_U_sets+1:end); %sd2x(2*n_U_sets)
sd2=size(D2,1);

%Flip D2
D2=[D2(:,n_U_sets+1:2*n_U_sets) D2(:,1:n_U_sets)]; %first n_U_sets-tuple is <= than the second

%Combine D1 and D2
D=[D1;D2];
sd=sd1+sd2; 

cr_increas=kron((1:1:sd), ones(1,n_vertices)); %1x[sd*n_vertices]
%n_vertices comes from the number of all possible vertices generated by each pair of n_U_sets-tuples

cc_increas=zeros(1,sd*n_vertices);
for j=1:sd
    D_temp=D(j,:);
    vertices=D_temp(C); %n_verticesxn_U_sets
    cc_increas(n_vertices*(j-1)+1:n_vertices*j)=changecoord_4(s1,s2, s3,idU1(vertices(:,1)), idU2(vertices(:,2)), idU3(vertices(:,3)), idU4(vertices(:,4))); %1xn_vertices
end

fill_increas=repmat(-C_sign, 1, sd);  %1x[sd*n_vertices]
%remember: signs are flipped because inequalities are <= 

RHS_increas=zeros(sd, 1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%% Compose %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Aeq=sparse([cr_obsequiv cr_zeros cr_one cr_symmetry], [cc_obsequiv cc_zeros cc_one cc_symmetry], [fill_obsequiv fill_zeros fill_one fill_symmetry], numeqconstraints, sU);   %[#equalityconstraints]x[sU]
Aineq=sparse(cr_increas, cc_increas, fill_increas, sd1+sd2, sU);   %[#inequalityconstraints]x[sU]
beq=[RHS_obsequiv; RHS_zeros;RHS_one; RHS_symmetry];


% MOSEK
prob.blx=zeros(sU,1); %lower bound unknowns
prob.ulx=ones(sU,1); %upper bound unknowns
prob.c=zeros(sU,1); %objective function
prob.a=[Aeq; Aineq];
prob.blc=[beq; -Inf(size(Aineq,1),1)];
prob.buc=[beq; RHS_increas];
[~,result]     = mosekopt('minimize echo(0)',prob, param_MOSEK);
end

